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1. 1  ntroduction 

1.1  Problem  Definition 

Energetic  systems  are  typically  designed  to  function  based  on  the 
occurrence  of  a  reaction  and  resulting  energy  release  in  a  body  of  energetic 
material.  Solid  explosives  are  typically  heterogeneous  bodies  composed  of  grains, 
voids  and  possibly  other  materials  such  as  binders.  A  useful  tool  in  studying 
energetic  systems  is  computational  modeling.  Modeling  allows  for  system  design 
while  minimizing  actual  manufacturing  and  testing.  In  the  case  of  weapon 
systems,  modeling  can  also  be  used  to  study  effectiveness  of  proposed  systems 
and  investigation  into  off-design  scenarios. 

Most  on-design  cases  for  weapon  systems  involve  a  detonation  occurring 
in  the  energetic  material.  In  terms  of  modeling  a  detonation,  there  are  two 
primary  ways  to  simulate  the  process.  The  first  is  to  prescribe  the  travel  and 
arrival  of  the  detonation  which  will  be  referred  to  here  as  a  prescription  method. 
The  second  is  to  solve  a  rate  law  equation  or  equations  to  simulate  the  reactions 
associated  with  the  detonation. 

1.2  Modeling  Options 

The  prescription  method  treats  the  shock  and  reaction  front  composing  a 
detonation  as  a  single  entity.  There  is  no  attempt  to  explicitly  represent  the 
reaction  process.  In  the  simplest  of  forms,  a  prescription  approach  maps  out  the 
arrival  of  the  detonation  wave  to  all  the  locations  throughout  the  energetic 
material  by  using  the  ignition  location  and  the  detonation  velocity  of  the 
material.  An  example  of  a  more  sophisticated  prescription  model  is  the 
Detonation  Shock  Dynamics  (DSD)  technique  that  has  been  used  to  successfully 
simulate  various  systems  ([3]).  It  is  typically  used  to  simulate  a  formed 
detonation  and  the  arrival  time  of  the  detonation  wave  to  the  various  parts  of 
the  energetic  material  are  mapped  out  using  methodologies  that  include  aspects 
such  as  the  curvature  of  the  wave.  It  has  been  applied  for  off-detonation  cases 
but  has  primarily  been  used  for  detonations  ([19]).  The  primary  limitation  of  the 
prescription  class  of  models  is  they  do  not  offer  the  universality  needed  to 
simulate  all  possible  occurrences  in  an  energetic  system. 

The  rate  law  approach  explicitly  models  the  reaction  process  within  the 
detonation  and  the  resulting  energy  release  which  drives  the  shock  in  the 
detonation.  There  are  a  variety  of  rate  law  type  models.  Some  are  similar  to  the 
multi-step,  finite-rate  temperature  based  Arrenhius  type  models  commonly  used 
in  gas  phase  combustion.  Others  employ  a  simple  one-step  model  with  the  rate 
being  dependent  on  the  local  pressure.  These  models  differ  in  their  formulation 
but  all  strive  to  represent  the  reaction  associated  with  the  detonation  as  a 
function  of  local  conditions  and  not  based  on  an  a  priori  mapping  of  the 
explosive  body.  In  this  study  various  rate  law  models  have  been  tested  and  each 
are  reviewed  in  the  following  sections  with  results  presented. 
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One  objective  of  the  current  study  is  to  determine  the  necessary  modeling 
methodology  to  capture  the  various  occurrences  that  are  of  interest.  The 
purpose  for  the  review  in  the  following  section  of  earlier  studies  is  to  determine 
what  sort  of  modeling  frame  work  is  needed  to  adequately  capture  all  the 
predominate  processes. 

1.3  Review  of  I  nitiation  and  Reaction  Mechanics 

When  dealing  with  detonation  modeling,  the  formulation  of  all  rate  law 
type  models  begin  with  an  investigation  into  the  processes  that  initiate  and 
sustain  the  reactions  ([23],  [17],  [28],  [8],  [16]).  It  is  universally  accepted  that 
the  reactions  result  from  "hot  spots"  that  appear  in  the  heterogeneous  reactive. 
The  genesis  of  these  areas  of  localized  energy  and  high  temperature  is  the  point 
of  debate  and  various  mechanisms  are  proposed  to  be  the  dominate  candidate. 
The  list  of  phenomena  "hot  spot"  formation  is  contributed  to  includes  jetting  of 
material  from  one  grain  to  another,  pore  collapse,  compression  of  gas  within  the 
voids,  friction  between  grains  and  plastic  deformation  of  the  material.  The 
relative  importance  of  the  various  mechanisms  that  can  cause  "hot  spot" 
formation  is  key  in  the  context  of  modeling  in  the  sense  of  the  mathematical 
formulation  of  the  model.  It  also  influences  the  resolution  needs  of  the  model 
which  affects  gridding  aspects  of  the  computational  framework. 

[28]  provides  a  short  review  of  the  various  mechanisms  to  which  "hot 
spot"  formation  is  attributed.  The  work  of  [15]  is  noted  to  conclude  that  gas 
compression  is  not  a  predominant  phenomena  for  hot  spot  creation  when  dealing 
with  shock  initiation.  The  role  of  viscoplastic  pore  collapse  in  driving  the  ignition 
phase  is  highlighted  as  is  the  role  of  burning  at  the  external  surface  of  the  grains 
in  the  growth  phase.  Modeling  is  performed  using  a  derivation  of  the  model  by 
[18].  [16]  uses  a  modeling  approach  that  does  not  explicitly  distinguish  the 
various  hot  spot  producing  mechanisms  but  takes  a  "unifying"  approach.  This 
assumes  the  hot  spots  can  be  treated  as  a  single  localized  heating  zone.  This 
approach  leverages  the  fact  that  in  the  end,  when  a  model  is  implemented  the 
finest  scale  at  which  differences  can  be  discerned  are  those  equivalent  to  the 
mesh  or  cell  space.  Therefore,  unless  resolution  is  achieved  at  scales  smaller 
than  the  grain  then  in  final  implementation  it  is  not  even  possible  to  distinguish 
the  various  hot  spot  production  mechanisms.  Therefore  a  "unified"  approach 
makes  sense.  The  current  study  follows  this  same  general  approach  and  details 
are  provided  later.  Here  a  brief  review  of  initiation  mechanisms  is  conducted  with 
the  objective  of  determining  if  the  modeling  framework  to  be  used  is  adequate  to 
represent  the  assumptions  on  which  the  various  rate  law  models  evaluated  are 
based. 

It  is  well  accepted  that  the  reaction  in  the  explosive  is  fundamentally 
associated  with  the  increase  in  temperature  in  the  material.  However,  when 
dealing  with  explosives  it  is  hard  to  isolate  the  temperature  rise  associated  with 
the  pre-ignition  phase  from  the  overall  temperature  increase  during  the 
detonation  event.  Therefore  a  review  of  the  response  of  non-reacting  materials 
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to  shocks  was  conducted.  Figure  1  shows  a  configuration  used  by  [4]  to  study 
the  effect  of  particle  shape  and  orientation  on  copper  compaction.  Both 
cylindrical  and  rectangular  particles  were  insulted  with  shocks  of  different 
strength.  Figure  2  shows  how  changing  the  shock  speed  and  therefore  strength 
affects  the  final  configuration  of  the  compacted  particles.  There  are  clearly 
noticeable  trends  as  the  shock  strength  is  increased.  Figure  3  shows  that  there 
are  material  jets  created  for  all  impact  cases.  These  jets  fill  the  voids  between 
the  particles  and  penetrate  the  particles  that  are  downstream  from  the  void.  The 
higher  velocities  result  in  deeper  penetration  and  more  deformation  of  all  the 
particles. 


Velocity  Boundary  Condition 


Transmitting  Boundary  Condition 

Figure  1.  Configuration  of  copper  particles  insulted  with  a  shock. 
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2.00  km/s 


Figure  2.  Effect  of  shock  strength  on  final  configuration  of  the 
compacted  particles. 


8 


Figure  3.  Noticeable  trends  as  shock  strength  is  increased. 

The  shock  strength  also  affects  the  temperature  rise  in  the  particles  which 
is  of  particular  interest  when  thinking  of  the  shock  compaction  process  in 
heterogeneous  explosives.  For  the  .25  km/s  case  the  temperature  of  the  interior 
of  the  particles  is  about  350K  while  the  exterior  is  380K  to  400K.  So  even  at  the 
lower  velocity  there  is  heating  of  the  surface  due  to  the  impact  of  the 
neighboring  particles  and  the  deformation  of  the  particle.  When  the  shock 
velocity  is  raised  to  1  km/s  there  is  a  considerably  different  temperature  field 
produced.  Figure  4  shows  the  variation  in  temperature  along  the  cut  line  denoted 
in  the  figure.  First  the  interior  of  the  particles  are  found  to  rise  to  about  500K 
and  the  temperature  at  the  surfaces  of  the  particles  reach  up  to  1500K  and 
greater. 
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1 .00  km/s 


Figure  4.  Distribution  of  temperature  across  the  particle  bed  for  the  1 
km/  s  case. 


Benson  [4]  also  investigated  the  effect  of  grain  shape.  The  configuration 
in  Figure  5  was  insulted  with  the  same  range  of  shock  strength  as  the  circular 
particles.  The  grain  shape  is  found  to  affect  the  load  path  of  the  stress  as  well  as 
the  jetting  of  the  material  as  evident  in  Figure  6.  The  range  of  temperatures 
resulting  from  the  compaction  process  are  the  same  as  the  cylindrical  particle 
cases.  Peaks  again  are  associated  with  material  jetting  and  plastic  flow  and 
aligned  particles  act  as  a  single  particle  with  heating  along  the  boundary  of  the 
grouping  as  evident  in  Figure  7.  Both  the  cylindrical  and  rectangular  particles 
show  that  deformations  and  heating  are  localized  at  the  surface  of  the  particles. 
Furthermore  the  heating  is  non-uniform. 


Figure  5.  Orientation  of  rectangular  particles  in  compaction  study. 
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Figure  6.  J  etting  of  material  and  stress  paths  in  the  case  of  rectangular 
particles. 
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1 .00  km/s 


Figure  7.  Heating  pattern  for  the  1  km/  s  case  for  rectangular  particles. 


Williamson  [39]  showed  that  energy  deposition  during  the  compaction  of 
particles  is  non-uniform.  For  circular  particles  undergoing  large  plastic 
deformation  the  frictional  heating  is  small.  And  due  to  the  short  compaction 
times  heat  conduction  and  radiation  are  negligible.  Using  the  configuration  of 
Figure  8  Williamson  [39]  determined  the  effect  of  shock  strength  on  the  final 
compaction  of  the  particles  and  associated  heating.  Figure  9  shows  the  impact 
velocity  affects  the  jetting  of  the  material  into  neighboring  particles  and 
influences  the  local  deformation.  Impact  regions  are  oriented  with  the  interstitial 
closures.  Figure  10  shows  the  change  in  heating  through  the  main  body  of  each 
particle  as  the  strength  of  the  shock  increases.  And  the  heating  is  found  to  be 
elevated  at  the  point  of  material  jetting.  The  temperature  time  histories  of  Figure 
11  show  the  large  spikes  in  temperature  occur  only  at  the  surface  of  the  particles 
and  correspond  to  the  void  closure. 


Figure  8.  Particle  configuration  used  to  evaluate  heating  during 
compaction. 
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Figure  9.  Final  particle  configuration  after  compaction  for  different 
strength  shocks. 
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Figure  10.  Heating  associated  with  the  compaction  process, 
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Figure  11.  Temperature  time  histories  at  a  location  internal  to  a 
particle  and  at  the  surface  of  a  particle. 
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The  role  of  the  void  is  explored  by  Williamson  as  well  [39].  Simulation 
results  such  as  that  in  Figure  12  show  the  clear  role  of  the  void  in  the  production 
of  hot  spots  in  a  shocked  heterogeneous  medium.  If  the  void  is  minimized  by  the 
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insertion  of  a  smaller  particle  the  temperature  rises  during  the  shocking  are 
greatly  reduced.  This  matches  the  data  from  [4]  that  shows  no  significant 
temperature  rise  in  bodies  that  are  in  contact  prior  to  the  arrival  of  the  shock. 
One  conclusion  from  these  earlier  studies  is  that  one  predominate  mechanism 
producing  hot  spots  is  the  impact  of  material  from  one  particle  (or  grain)  with 
another  particle  (or  grain). 
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Figure  12.  Effect  of  particle  configuration  in  the  compaction  and 
resulting  heating. 
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Figure  13.  Effect  of  pore  gas  in  the  resulting  heating  of  compacted 
particles. 


As  mentioned,  the  role  of  the  compression  of  the  gas  found  in  the  pores 
in  driving  the  reactions  is  routinely  cited.  Some  insight  into  the  aspect  of  the 
problem  can  be  shed  by  considering  the  early  work  by  [39]  where  the  heating  of 
heterogeneous  materials  were  simulated  assuming  an  air  filled  pore  as  well  as  an 
evacuated  pore.  Figure  13  shows  when  the  pore  is  evacuated,  the  heating  of  the 
particle  follows  a  certain  pattern  due  to  material  compressions  and  resulting 
heating.  When  the  pore  is  filled  with  air,  the  pattern  is  the  same  but  the  level  of 
heating  is  increased.  This  is  due  to  the  fact  that  with  air,  and  for  the  same  travel 
of  the  compressing  plate,  the  grain  must  be  compressed  more  and  therefore  the 
deformation  is  higher  resulting  in  higher  temperatures.  Therefore  the  air  is  acting 
as  additional  mass  impacting  the  grain.  For  the  case  of  energetic  material,  this 
heating  can  initiate  reactions. 

There  is  a  difference  between  the  compaction  of  the  inert  materials 
reviewed  here  versus  energetic  materials.  In  the  simulations  by  [39]  there  is  not 
the  potential  of  reactions  occurring  between  gasified  particle  material  and  the 
pore  air.  This  potential  process  would  be  accurately  represented  by  capturing  the 
melting  at  the  grain  surface  and  resulting  gasified  fuel  mixing  with  pore  air  and 
resulting  gas  phase  combustion. 

The  scenario  studied  by  Williamson  [39]  was  reproduced  using  the  code 
cutter  provided  to  us  by  the  Munitions  Directorate.  In  the  simulation  a  1-D 
configuration  was  considered  and  three  grains  of  copper  were  used.  The  material 
model  for  the  grains  was  the  Mie-Gruenison  model.  The  grains  were  compressed 
as  a  rate  of  1  km/s  and  2  km/s.  The  resulting  heating  patterns  are  shown  in 
Figure  14.  The  trends  are  similar  to  those  seen  in  the  earlier  work  by  [39]  with 
heating  being  concentrated  at  the  surface  of  the  grains.  The  magnitudes  are 
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similar  as  well. 

Simulations  using  cutter  were  also  run  for  three  grains  of  explosive  to 
simply  test  the  modeling  framework.  Properties  for  TNT  were  used  but  no 
reactions  were  allowed  to  occur.  The  predicted  compaction  and  resulting  heating 
patterns  early  and  late  in  the  event  can  be  seen  in  Figure  15.  The  results  clearly 
show  early  in  the  event,  the  compaction  at  the  surface  of  the  grains  produced 
localized  heating  as  seen  in  the  inert  materials.  However,  the  magnitude  of  the 
heating  is  rather  low.  Also  shown  in  the  figure  is  the  heating  history  of  the  grains 
with  each  line  denoting  a  different  point  in  time.  The  heating  is  found  to  not 
remain  concentrated  at  the  grain  boundaries  but  is  found  to  be  elevated  at  the 
interior  of  the  grain.  This  does  not  follow  the  expected  behavior  for  an  "inert" 
material  and  is  most  likely  due  to  the  inadequacy  of  the  constitutive  model  used. 
More  on  the  exact  modeling  used  in  this  study  is  provided  later. 

It  is  known  that  various  mechanisms  can  lead  to  the  initiation  of  reactions 
in  heterogeneous  explosive  materials  such  as  void  collapse  and  material  jetting. 
It  is  clear  based  on  the  review  of  micro-scale  studies  involving  inert  materials 
that  the  hot-spot  formation  is  concentrated  at  the  surface  of  the  grains.  And  in 
the  configurations  studies,  the  primary  heating  comes  from  material  of  one  grain 
impacting  neighboring  grains.  The  impact  causes  deformation  of  the  grain  and 
associated  heating.  Earlier  studies  have  supported  the  conclusion  that  plastic 
deformation  is  the  dominate  feature  of  the  energy  localization  [2].  Two 
observations  can  be  made  from  these  earlier  studies.  First,  to  advance 
alternative  modeling  options  for  explosives,  appropriate  constitutive  models  are 
needed  for  explosive  materials.  The  heating  of  the  grain  due  to  material 
deformation  needs  to  be  accurately  resolved.  The  data  also  suggests  that  the 
dynamics  at  the  smaller  scales,  i.e.  grain-to-grain  impact  and  resulting 
deformation,  may  be  adequately  represented  without  exactly  resolving  down  to 
the  individual  grain.  Since  heating  is  focused  at  the  grain  surface,  a  modeling 
framework  that  accounts  for  parameters  such  as  volume  fraction  of  solid  and 
gaseous  material  and  estimates  the  surface  area  may  be  adequate  to  perform 
simulations. 


Figure  14.  The  compaction  and  resulting  heating  of  copper  particles  as 
predicted  by  cutter. 
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Following  the  initiation  of  burning  in  an  energetic  material  is  the  growth  of 
the  reaction.  The  growth  of  reactions  in  any  material  are  tied  to  the  material 
properties.  It  is  also  tied  to  the  size  and  shape  of  the  burning  front.  For  the 
scenarios  of  interest  here  the  properties  are  of  course  those  of  the  solid  material 
and  the  gases  produced  from  the  solids.  The  size  and  shape  of  the  burning 
fronts  will  be  associated  with  the  grain  sizes  and  orientation  since  the  burning 
originates  on  the  surface  of  the  grain.  The  evolution  of  the  fronts  will  depend  on 
the  post-compaction  configuration  of  the  grains  and  properties  of  the  material 
such  as  the  ablation  properties  of  the  grains.  It  is  obvious  that  a  micro-scale 
model  with  adequate  resolution,  material  properties,  and  chemical  reaction 
models  would  accurately  predict  the  growth  of  the  reactions.  The  question  is  can 
a  meso-scale  model  better  represent  the  micro-scale  processes  than  some  of  the 
macro-scale  models  currently  used. 


SwniAHion  THT  grain  impaction 
Pane*?  Mfltwny  t  Kmte 
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Figure  15.  The  compaction  and  resulting  heating  of  TNT  particles  as 
predicted  by  cutter. 
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2.  Rate  Law  Models 

Here  we  will  refer  to  any  reaction  model  that  uses  local  parameters  such 
as  pressure,  temperature  and  material  concentrations  to  specify  a  reaction  or 
burning  rate  as  a  rate  law  model.  Figure  16  shows  the  relative  relationship 
between  grains  and  cell  size  when  using  the  various  models  reviewed  here. 


Figure  16.  Configuration  of  multiple  grains  within  a  computational  cell 
or  element. 

2.1  Ignition  and  Growth 

One  popular  rate  law  model  is  the  ignition  and  growth  model  by  [23]. 
The  model  uses  mixed  elements  or  volumes  that  contain  unreacted  explosive  and 
reaction  products  at  equilibrium  in  terms  of  pressure.  The  temperature  of  the 
two  materials  are  not  in  equilibrium  and  the  difference  in  whether  temperature 
or  pressure  equilibrium  is  assumed  was  quantified  in  [11].  The  model  framework 
uses  the  J  ones-Wilkins-Lee  (JWL)  equation  of  state  for  both  materials.  The 
formulation  assumes  a  small  fraction  of  the  explosive  is  ignited  by  the  passing 
shock  in  the  detonation.  This  follows  the  experimental  evidence  of  the  localized 
hot-spot  formation.  The  growth  is  modeled  to  be  driven  by  the  pressure  build  up 
in  the  material  and  the  surface  area  of  the  burning  front.  Using  the  assumption 
of  micron-sized  spherical  burning  regions  that  grow  outward,  the  rate  of  energy 
release  is  modeled  using 


dFdt  =  /(l  -  F)x  +  G(l  -  F)x  FyPz 
V  =  V0!Vl-\ 


(1) 


where  F  is  the  fraction  of  the  explosive  that  has  reacted  and  t  is  time.  The  initial 
specific  volume  of  the  explosive  is  denoted  by  V0  and  that  of  the  shocked  but 

unreacted  explosive  by  vx .  The  pressure  is  denoted  by  P  (in  megabars)  and  I,  x, 
r,  G,  y,  and  z  are  constants.  The  various  terms  in  the  equation  and  associated 
constants  are  used  to  mimic  aspects  such  as  the  coupling  between  compression 
of  material,  pressure  and  amount  of  plastic  work  involved.  The  constant  G  is 
related  to  the  surface  area  to  volume  ratio  for  the  particular  material  and  the  Pz 
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incorporates  a  pressure  dependent  laminar  burning  rate.  The  term  (l-F)xFy 
ensures  the  rate  zeros  when  there  is  no  more  unreacted  material  and  the 
exponents  set  the  value  of  F  at  which  the  growth  of  the  reaction  is  maximized. 

Results  using  this  model  will  be  presented  here  which  show  the  model 
correctly  captures  the  dynamics  of  a  detonation.  However,  earlier  work  suggests 
the  model  has  limited  applicability  in  non-detonation  cases  [12]. 

2.2  Hot  Spot  and  Balance 

Another  rate  law  model  was  proposed  by  [17]  where  the  combustion 
process  is  explicitly  divided  between  the  hot-spot  region  and  the  rest  of  (or  the 
balance  of)  the  energetic  material.  The  model  uses  temperature  dependent 
induction  times  and  explicitly  tracks  intermediate  and  final  material  states.  A 
continuum  approach  is  used  to  represent  each  element  or  cell.  Studies  used  to 
formulate  the  model  in  part  by  studying  the  fact  that  the  temperature  found  in 
hot-spots  correlates  well  with  measured  induction  times  in  explosives  such  as 
HMX.  The  model  framework  uses  a  mixed  cell  type  approach  with  each 
component  at  its  particular  state.  So  within  a  given  cell,  there  is  a  certain  fraction 
of  material  that  is  at  the  hot-spot  state  while  some  is  at  an  intermediate  state 
and  the  remainder  is  at  the  final,  post  reaction  state.  If  the  assumption  is  made 
that  within  a  cell,  all  hot-spots  are  at  the  same  temperature  simulation  results 
capture  the  peak  values  in  the  time  history  of  properties  such  as  material 
velocities  but  the  actual  time  history  is  not  correctly  resolved.  In  this  mode,  the 
intermediate  states  are  essentially  ignored.  If  the  intermediate  states  are 
introduced  then  the  simulations  are  found  to  correctly  capture  the  time  history  of 
the  detonation  process. 

The  hot-spot  and  balance  approach  is  generally  the  same  approach  as  the 
ignition  and  growth  in  that  the  different  states  of  the  material  are  represented 
from  unreacted  solid  energetic  material  to  gaseous  detonation  products.  It  is 
different  in  that  the  rate  of  energy  release  is  temperature  based  with  the  form  of 
the  rate  equation  being 


f  =  (l-f)Zexp[~a(0+fif)]  (2) 

where  a ,  (5 ,  and  Z  are  constants  for  the  explosive  of  interest  and  6  is  the  local 
temperature.  In  the  case  of  both  models,  empirical  parameters  are  needed  to 
close  the  models  and  the  models  are  measured  and  applied  at  a  continuum  level 
that  incorporates  both  the  material  grains  and  pores. 

2.3  Temperature  Based  Kinetics 

The  models  briefly  reviewed  in  the  previous  two  sections  were  formulated 
with  the  main  objective  being  to  simulate  detonations  using  local  conditions  as 
opposed  to  using  a  prescribed  type  model.  The  detonation  scenario  represents 
one  end  of  the  spectrum  in  terms  of  time  scales  of  mechanical  and  thermal 
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processes.  At  the  other  end  of  the  spectrum  in  terms  of  time  scales  is  the 
problem  of  thermal  decomposition  of  explosive  materials.  The  problem  of 
thermal  decomposition  during  scenarios  such  as  cook-off  are  routinely  modeled 
using  multi-step  kinetics  models  where  intermediate  products  are  tracked  and 
the  rates  are  dependent  on  temperature.  It  is  a  goal  to  eventually  model  the 
impact  and  shock-induced  reactions  using  the  same  modeling  approach  [36]. 
This  option  will  be  explored  here  using  a  model  proposed  in  [40].  In  this  earlier 
study,  scaled  thermal  explosion  experiment  (STEX)  scenarios  are  modeled  using 
a  multi-step,  multi-species  reaction  mechanism. 

In  the  work  by  [40],  a  three-step,  four-species  mechanism  for  C-4  of  the 

form 


B,r |  =  Zpxp{-  EJRT )pA 

(3) 

B  —>  C,  r2  =  Z2exp{-  EJRT )pB 

(4) 

C  — »  D,r3  =  Z:exp(-  EJRT  )p^, 

(5) 

where  A  and  B  are  solid  species,  and  C  and  D  are  gases.  The  parameter  rt  is  the 
reaction  rate  and  Z;  is  the  rate  factor  and  Et  is  the  activation  energy  for 

reaction  i .  This  form  is  similar  to  the  Arrenhius  rate  law  typically  used  in  gas 
phase  combustion  modeling.  The  parameter  pj  is  the  mass  fraction  of  species 

j  ■ 
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3.  Geometry  Representation 

There  exists  a  variety  of  meshing  frameworks  for  use  in  computational 
modeling.  The  main  classes  are  (1)  body-fitted  curvilinear  meshing,  (2) 
unstructured  tetrahedral  meshing,  and  (3)  Cartesian  meshing.  All  have  their  pros 
and  cons. 

Because  the  problems  of  interest  here  can  range  in  size  and  configuration 
it  was  decided  to  use  the  Cartesian  meshing  approach  with  the  addition  of  mesh 
adaptation  and  refinement.  The  method  employed  here  uses  the  basic  principles 
found  in  the  work  by  [7],  [5],  and  [6].  A  short  review  is  provided  here  but  the 
reader  is  referred  to  these  earlier  works  for  more  details. 

Using  the  Cartesian  approach,  the  domain  of  interest  is  defined  as  any 
rectangular  shape.  Also  defined  is  the  size  of  the  coarsest  grid  desired  which  is 
referred  to  as  the  Level  0  cells.  The  cells  need  not  be  cubic  but  larger  aspect 
ratio  cells  can  introduce  numerical  error  depending  on  the  particular  flow  solver 
used.  Also  defined  is  the  number  of  adaptation  levels  desired,  referred  to  here  as 
njevels.  The  two  parameters,  Level  0  cell  size  and  njevels  will  influence  both 
solution  accuracy  and  efficiency.  This  will  depend  on  the  particular  problems 
being  solved  and  benchmark  cases  will  be  solved  here  to  determine  appropriate 
settings  for  the  problems  of  interest. 

The  Cartesian  mesh  used  here  will  be  adapted  as  needed  to  define  either 
geometries  or  solutions.  The  particulars  of  the  adaptation  will  be  presented  in 
following  sections.  Here  the  fundamentals  of  the  meshing  structure  are 
discussed.  The  domain  of  interest  is  decomposed  using  the  specified  Level  0  cells 
as  depicted  in  Figure  17.  If  increased  resolution  is  needed  these  cells  are 
systematically  divided  into  8  children.  The  requirement  for  refinement  is 
dependent  on  a  local  feature  of  interest  in  the  geometry  or  solution.  In  the 
current  implementation,  the  discontinuity  between  neighboring  cells  is  never 
greater  than  a  2:1  ratio.  This  restriction  does  require  that  some  additional  cells 
be  refined  other  than  just  those  meeting  specified  criteria.  However,  it  simplifies 
the  flow  solver  and  difference  equation  solution  algorithms  and  more  than 
adequately  accounts  for  any  additional  cost. 


21 


*x 


Figure  17.  The  cell  division  used  when  refinement  is  required.  A  cell  at 
level  /7is  divided  into  8  children  at  level  n+1. 

To  facilitate  the  management  of  the  computational  mesh  and  the  solution 
process,  there  are  some  basic  grid  accounting  parameters  that  need  to  be 
associated  with  each  cell.  These  include  the  following: 


•  level  =  level  of  the  cell 

•  idiv=  0  if  cell  is  not  divided,  =  1  if  cell  is  divided 

•  Parent  =  pointer  to  a  parent  cell 

•  C2[i][j][k]  =  pointer  to  children  cells:  i=0,l;  j =0,1;  k=0,l 

•  xN  =  pointer  to  a  neighboring  cell:  ;r=L(eft),  R(ight),  B(ottom),  T(op), 
U(nder),  O(ver) 


The  first  three  parameters  are  self  evident.  The  third  parameter  defines 
the  children  cells  if  any  cell  at  any  level  is  divided.  The  pointers  are  stored  in  a 
three-dimensional  array  making  it  easy  to  define  algorithms  such  as  flow  solvers. 
Figure  18  shows  the  definition  of  the  cells  and  orientation  of  the  numbering.  The 
symbol  #  denotes  that  from  that  vantage  point  the  index  could  be  0  or  1. 
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Figure  18.  Definition  of  the  children  cells. 


4.  Computational  Model 


Here  the  computational  model  used  to  simulate  the  scenarios  of  interest  is 
detailed.  The  wave  propagation  and  material  transport  processes  found  in  the 
scenarios  of  interest  can  be  modeled  using  the  Euler  equations  which  can  be 
written  in  conservative  form  as 


M+M(eLn(e) 

dt  dx,  ; 


(6) 


where  the  subscript  /'denotes  direction  and  the  vectors  are 
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(7) 


where  the  subscript  /  is  1,  2  or  3  respectively  to  denote  the  directions  x,  y  or  z. 
The  operator  Sm  works  such  that  the  pressure  flux  functions  only  on  that 
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component  of  momentum  coinciding  with  the  direction  of  the  flux.  The 
parameters  p ,  u,  v,  and  w  are  respectively  the  density,  x,  y,  and  z  components 
of  velocity  and  P  is  pressure.  The  mass  fraction  of  material  n  is  denoted  by  an 
and  con  is  the  source  term  for  that  particular  material.  The  parameter  NS  denotes 

the  number  of  different  materials  (or  species)  involved  in  the  problem  and  since 
the  overall  mass  (or  density)  is  conserved  then  the  mass  fraction  of  only  NS-1 
species  have  to  be  explicitly  tracked  since  the  rules 

NS  NS 

P  =  YfinP^  IX  =1  (8) 

n— 1  n— 1 


hold.  The  source  term  for  each  material  will  be  discussed  in  detail  in  a  following 
section. 

The  parameter  H  represents  the  total  enthalpy  (H  =  E  +  P/p)  and  total 
energy  is  represented  by  the  parameter  E  and  is  the  sum  of  internal  (e)  and 

kinetic  (-X2)  energy.  For  the  multi-material  problem,  the  internal  energy  is  a 
sum  of  contributions  from  each  species  through  the  equation 


e  = 


NS 

Yjanen 


n= 1 


(9) 


and  the  representation  of  Pn  depends  on  the  nature  of  material  n.  For  the  case 
of  an  ideal  gas  the  partial  pressure  can  be  represented  using 

P„=pR,T/M„  (10) 

n 

where  Ru  is  the  universal  gas  constant,  T  is  temperature  and  Mw  is  the 

wn 

molecular  weight  of  material  n. 

For  the  energetics  problems  of  interest  here,  materials  involved  such  as 
the  solid  explosive  and  detonation  products,  require  a  non-ideal  equation  of  state 
to  be  used.  Here  the  Jones-Wilkins-Lee  (J  WL)  equation  of  state  will  be  used  and 
can  be  found  in  two  forms,  an  internal  energy  dependent  form 

P  =  A(  1  -  colRxV)eRxV  +  5(1  -  colR2V)eR lV  +  cupel  V  (11) 

and  a  temperature  dependent  form 


-rV  -rV 

P  =  Ae  1  +Be  2  +a>CvT/V. 


(12) 
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The  parameters  A,  B,  /?, ,  i?2  and  are  constants  associated  with  the 
material,  Cv  is  the  average  heat  capacity  and  V  is  the  relative  volume,  pjp , 
with  pa  being  the  initial  reference  density. 

When  solving  a  problem  involving  multiple  materials,  the  set  of  equations 
must  be  closed  assuming  either  the  materials  are  in  equilibrium  in  terms  of 
pressure  or  temperature.  For  the  high  explosive  scenarios,  the  standard 
assumption  is  pressure  equilibrium  since  pressure  can  equalize  much  faster  than 
temperature.  Using  the  mass  fraction  approach  and  the  equilibrium  pressure 
assumption,  the  contribution  from  the  non-ideal  materials  to  total  internal  energy 
can  be  represented  using 


where 


pei  =  VPIa>i+Ai+Bi-E0' 


4  =  4 


B,  =  B, 
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R,  o.)i 
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e  1 


(13) 


(14) 


The  energy  component,  Ea,  is  analogous  to  the  heat  of  formation  of  the 

material  and  since  the  convention  commonly  used  with  the  J  WL  is  to  define  this 
quantity  as  a  positive  number,  the  value  must  be  subtracted  from  Equation  (13) 
to  produce  the  correct  energy  release  effect.  The  pressure  of  the  mixture  can  be 
calculated  using 


p  = 


NS  /  ^  \ 

/*J-&l4  +  4-4j 

i= 1 


NS 


i= 1 


-1 


and  the  temperature  of  each  material  can  be  calculated  using 

T  =^cr[p-Aiexp(-Rlv)-Biexp(-R1v\ 


(15) 


(16) 


The  temperature  of  the  mixture  is  then  calculated  using  a  mass  weighted 
average 


NS 


t  =  Yp?, 


1=1 


(17) 


25 


The  values  of  the  constants  are  presented  later  along  with  the  results. 

There  are  a  variety  of  ways  to  represent  the  flux  terms  F  and  two 
alternative  methods  were  evaluated  here.  Details  of  these  methods  will  be 
presented  in  a  following  section.  The  solution  of  the  source  terms  is  also 
discussed  in  a  following  section. 

When  considering  the  solution  of  the  set  of  equations  and  the  integration 
of  the  equations  through  time,  there  are  also  a  variety  of  options  with  the  two 
main  categories  being  implicit  and  explicit  methods.  For  the  high-speed  problems 
of  interest  here,  the  stability  advantages  of  implicit  methods  which  allow  for 
large  time  steps,  are  not  realized  in  that  small  time  steps  are  needed  to  resolve 
the  phenomena  of  interest  [10].  However,  the  reaction  source  term  introduces  a 
stiff  equation  and  using  implicit  methods  do  have  an  advantage.  Therefore,  here 
a  splitting  strategy  called  Strang  splitting  that  has  been  used  successfully  in  high¬ 
speed  reacting  flows  is  used.  The  dependent  variables  Q  are  advanced  from 
state  n  to  n+1  by  successively  applying  the  fluid  dynamic  operator  (lf)  and 
reaction  operator  ( Ln )  in  the  order 


Qn+l  =LAfL*LAfQ"  (18) 

The  nature  of  the  fluid  dynamics  and  reaction  operators  depend  on  the 
characteristics  of  the  models  being  used.  For  the  work  here,  the  fluid  dynamics 
operator  used  will  be  explicit  while  the  reaction  operator  will  be  explicit  or 
implicit  depending  on  the  specifics  of  the  reaction  model  being  used. 
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5.  Numerical  Comparisons 


5.1  2D  Scenario  Comparisons 

Because  the  objective  here  is  to  move  towards  simulating  energetic 
systems,  the  two  different  numerical  flux  schemes  were  also  compared  within 
the  framework  of  the  adaptive  meshing  to  see  if  either  method  demonstrated 
any  noticeable  advantages.  The  test  scenario  was  a  simplified,  post  detonation 
state  mimicking  the  flyer  plate  scenario  shown  in  Figure  92.  A  volume  equivalent 
to  the  explosive  was  set  to  the  post  detonation  state  in  terms  of  material, 
pressure  and  density  as  depicted  in  Figure  93.  In  the  figure  are  pressure  (P), 
density  (r),  x-velocity  (u),  y-velocity  (v)  and  mass  fraction  of  the  detonation 
product  (DP)  plotted  along  with  the  dynamic,  adapting  grid.  Also  plotted  is  the 
adaptation  criteria  (s)  and  the  level  to  which  the  local  grid  is  resolved  (level). 

Figure  94  shows  the  results  using  the  two  different  methods  at  a  point  in 
time  shortly  after  the  initiation  of  the  simulation  while  Figure  95  compares  the 
results  some  time  later.  Both  schemes  performed  well  within  the  adaptive  mesh 
frame  work.  There  were  also  no  significant  differences  in  the  results.  There  is 
some  difference  in  the  contour  coloring  in  the  pressure  and  density  plots  but  this 
is  attributed  to  the  scaling  during  the  plotting  and  not  due  to  differences  in  the 
results. 


mylar  flyer  plate 


Explosive  (1 8mm  x  9mm)  target  plate 


Figure  92.  Impact  scenario  used  to  compare  the  numerical  flux 
methods. 
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Figure  94.  Comparison  between  the  Steger- Warming  (left)  and  AUSM 
(right)  results  shortly  after  initiation  of  the  simulation. 
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Figure  95.  Comparison  between  the  Steger- Warming  (left)  and  AUSM 
(right)  results  some  time  later  after  initiation  of  the  simulation. 

5.2  Findings 

In  both  the  ID  and  2D  simulations,  there  are  no  significant  differences 
found  between  the  results  from  the  Steger-Warming  scheme  and  the  AUSM  and 
AUSM+  schemes.  There  are  some  boundary  issues  associated  with  AUSM  but 
those  can  be  resolved  using  the  approaches  mentioned.  Because  the  AUSM 
schemes  do  not  represent  any  noticeable  advantages  we  proceeded  with  using 
the  Steger-Warming  scheme  because  it  has  demonstrated  excellence  robustness 
in  previous  work. 
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6.  Reaction  Source  Term  Representation 

This  section  contributed  by  Andrew  Staszak  and  J .  Keith  Clutter 

The  governing  equations  have  been  presented  and  the  reaction  source 
term  is  critical  to  representing  the  explosion  processes.  The  options  in  how  to 
model  these  processes  were  outlined  earlier.  All  these  models  constitute  a  set  of 
equations  that  have  to  be  solved.  Here  we  evaluate  different  options  to  solve  the 
equations.  The  most  complicated  model  considered  here  for  the  energetic 
materials  problem  is  the  multi-step,  finite-rate.  The  reaction  models  in  these 
systems  are  very  complex  and  require  a  great  deal  of  computational  effort.  When 
considering  chemical  reaction  simulations,  two  main  factors  need  to  be 
considered:  the  speed  of  the  reaction  and  the  characteristic  time  scale  of  the 
flow  itself.  In  applications  where  the  flow  field  is  fast  moving  a  multi-step,  finite- 
rate  reaction  model  must  be  used  to  model  the  combustion  process. 

Current  methods  used  to  solve  combustion  equation  systems  have  many 
issues  with  the  extensive  computational  time  needed  to  compute  the  chemical 
kinetics.  This  computational  complexity  arises  because  many  of  the  systems 
characterized  in  these  types  of  problems  are  characterized  on  a  wide  range  of 
temporal  and  spatial  scales.  This  range  difference,  in  most  cases,  is  many  orders 
of  magnitude.  This  range  of  scales  is  manifested  as  "stiffness"  in  the  sets  of 
differential  equations.  Because  of  this,  many  industry  models  have  been  reduced 
to  using  complex  combustion  systems  only  in  simple  flow  problems  [9].  To 
counter  these  inherent  extensive  calculations,  much  of  the  current  effort  is 
centered  around  mechanism  reduction  methods,  which  reduce  the  number  of 
reactions  in  the  mechanism  to  a  few  primary  steps.  To  reduce  mechanisms  in 
this  manner  requires  investigative  time  into  the  behavior  of  each  reaction  in  the 
full  mechanism.  Reduction  methods  have  been  noted  to  have  three  major 
drawbacks  [27]: 

1.  Reductions  are  specific  to  a  mechanism  and  require  additional  human 
time  and  labor  to  develop. 

2.  Partial-equilibrium  and  steady-state  assumptions  can  only  be  applied  to 
specific  ranges  in  composition,  pressure,  temperature,  and  time  scales.  Outside 
of  these  ranges  results  can  be  poor. 

3.  Error  estimates  can  only  be  obtained  by  comparing  the  solution  of  the 
detailed  mechanism  to  that  of  the  reduced  mechanism. 

Alternative  methods  which  can  reduce  computational  intensity  in  existing 
finite-rate  chemistry  models  would  allow  for  greater  ease  in  the  modeling  of 
complex  combustion  mechanisms  in  high  speed,  complex,  reacting  flows.  In  this 
paper,  the  implementation  of  numerical  derivatives  in  solving  the  governing 
equations  is  evaluated.  The  proposed  approach  is  compared  to  solution  methods 
commonly  used  in  the  simulation  tools  for  reactive  flows.  The  three  hydrogen- 
oxygen  reaction  mechanisms  shown  in  Figure  96,  97,  and  98  will  be  used  to 
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perform  the  comparisons.  This  will  provide  a  range  of  differing  complexity  to 
capture  any  performance  dependence  on  mechanism  size. 
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Figure  96.  Seven-step  Hydrogen- Oxygen  Reaction  Mechanism. 
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Figure  97.  Nineteen-step  Hydrogen- Oxygen  Reaction  Mechanism. 
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Figure  98.  Thirty  two  -  step  Hydrogen- Oxygen  Reaction  Mechanism. 

6.1  Governing  Equations 

The  governing  equations  for  the  complete  propulsion  system  problem 
include  the  Navier-Stokes  equations.  Because  the  focus  here  is  on  the  chemical 
reaction  aspects,  the  governing  equation  can  be  represented  in  the  vector  form 

=  Q  (19) 

ot 

where  only  the  dependent  variables  and  their  source  terms  remain  in  the 
equation.  Equation  (19)  represents  the  constant  density  chemical  reaction  model 
that  will  be  used  in  this  paper.  The  terms  Q  and  Q  represent  vectors  containing 
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density  and  terms  specific  to  each  species  represented  in  a  mechanism.  The 
vector  of  dependent  variables  is 


Q  =  [p,pal,---,paNS_J  (20) 

where  a  is  the  mass  fraction  of  the  ith  species,  and  NS  is  the  total  number  of 
species  in  the  mechanism.  The  source  terms  are  represented  in  the  vector 

Q  =  (21) 


where  the  terms  are  only  represented  to  the  NS-1  species  since  overall  mass 
conservation  is  maintained.  The  source  terms  represent  the  production  rates  of 
each  individual  species  in  the  mechanism.  These  terms  incorporate  both  the 
forward  and  backward  production  rates  of  each  species  throughout  all  of  the 
reactions  of  the  mechanism.  The  solution  assumes  a  mixture  of  perfect  gases 
and  conserves  total  energy.  The  specific  heat  of  each  species  is  represented 
using  J  ANAF  data  [29]. 

The  chemical  model  is  represented  by  a  reaction  mechanism  containing 
any  number  of  reactions  (NR)  and  species  (NS).  These  reactions  are  in  the  form 


NS 


NS 


Ir'.,x,kf  I/Ixj  ■  '  1 . v'( 


(22) 


J=l  k  J=1 

Kbi 


where  i /.  and  v.'  are  the  stoichiometric  coefficients  for  species  j  in  reaction  i, 

and  Xj  is  the  molecular  concentration  of  species  j.  The  forward  and  backward 
reaction  rates  are  represented  in  the  terms  kfi  and  kbi  for  each  reaction.  The 

forward  reaction  rates  are  calculated  using  the  Arrhenius  rate  expression 


kn  =  AT  ' exp 


RT 


(23) 


where  Ru  is  the  Universal  Gas  Constant.  Ai  is  the  pre-exponential  factor,  bi  is  the 

temperature  exponent,  and  Eai  is  the  activation  energy  for  the  forward  reaction 
and  values  used  can  be  found  in  Figures  96,  97,  and  98.  The  backward  rates  are 
calculated  using  the  relationship  between  the  equilibrium  constant  and  the 
forward  and  backward  rates.  The  production  rates  of  each  individual  species  are 
calculated  through  a  summation  of  contributions  from  each  reaction  and 
represented  symbolically  as 


m 

d>j  =  mwj^jy'.j  -v'y)- 


1=1 


NS 


NS 


kA K" 


n— 1 


n= 1 


(24) 
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Where  mwj  represents  the  molecular  weight  of  the  species. 

The  governing  equations  result  in  a  system  of  equations  that  must  be 
integrated  in  time.  When  coupled  with  the  fluid  dynamic  terms  the  time  step  is 
typically  set  by  the  CFL  condition  associated  with  the  flow.  The  restrictions  from 
the  fluid  aspects  of  the  problem  are  typically  less  restrictive  than  the  reaction 
aspects.  A  common  approach  is  the  use  of  a  split  operator  method  that  allows 
different  solutions  methods  to  be  used  for  the  fluid  dynamic  operation  versus  the 
reaction  operations  [10].  This  allows  the  overall  time  step  to  be  set  based  on 
parameters  of  the  study.  If  it  is  set  such  that  the  restrictions  associated  with  the 
fluid  dynamics  will  define  the  time  step  between  state  n  and  n+1,  then  the 
reaction  operator  is  used  to  advance  the  chemical  reactions  from  one  state  to 
the  next.  The  focus  here  is  on  numerical  methods  used  for  the  reaction  operator. 
The  reaction  system  of  equations  can  be  solved  in  a  variety  of  means  and  one  of 
the  simplest  options  is  the  use  of  explicit  methods.  Explicit  schemes  allow  for  the 
dependent  variables  to  be  directly  computed  in  terms  of  known  quantities.  In  the 
reaction  model,  the  partial  derivative  of  the  dependent  variables,  Q,  with  respect 
to  time  can  be  found  using  a  first  order  upwind  difference  operator  [21].  In 
difference  form,  the  system  of  equations  is 

Qn+l=Qn+QnAt  (25) 

where  the  superscript  denotes  the  time  level  at  which  the  vector  of  variables  is 
represented.  The  At  is  chosen  based  on  the  reaction  mechanism  that  is  being 
modeled.  Although  explicit  schemes  are  simple,  they  have  very  low  stability  and 
can  easily  produce  non-physical  solutions.  Because  of  this,  careful  attention  must 
be  paid  when  choosing  a  step  size. 

Another  option  is  the  use  of  a  Runge-Kutta  scheme.  This  approach  is  an 
attempt  to  benefit  from  the  simplicity  of  explicit  methods  while  decreasing 
computational  time.  Here  the  Fourth-Order  Runge-Kutta  will  be  used.  The 
reaction  model  is  represented  as  a  first  order  differential  with  initial  values 

=  Q.(t,  pa),  pa(t0)=  pa°  (26) 

ot 

where  the  initial  values  would  be  Qn .  By  formal  integration  of  the  derivative  and 
applying  Euler's  algorithm,  the  fourth-order  method  can  be  written  as  [32] 
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K  =hf(tn,pan) 

K  =  hf^tn  +—,pan  +^- 

h=hf[tn+^,pccn  +  1^  (27) 

k4  =  hf[tn  +h,pan  +kj) 

pan+x  =  pan  +^  +  ^  +  ^  +  ^  +  o(A5) 

6  3  3  6 

This  method  requires  the  evaluation  of  the  source  terms  in  the  governing 
equations  four  times  for  each  time  step.  Thus,  the  value  at  the  next  time  step  is 
found  using  the  present  value  plus  a  weighted  average  of  the  slopes  at  several 
points.  The  goal  of  the  Fourth-Order  Runge-Kutta  method  is  to  increase  the  size 
of  time  step  that  can  be  taken  compared  to  the  generic  explicit  method. 

Another  option  is  to  solve  the  set  of  equations  using  an  implicit  approach. 
Such  an  approach  helps  address  issues  such  as  stability  in  the  solution  process 
and  numerical  restrictions  associated  with  the  time  steps.  The  difference 
equations  that  are  actually  solved  come  from  the  linearization  of  the  chemical 
reaction  source  term 


where 


Q"+1  =  nn+(A)nAQ 


(28) 


(29) 


is  the  Jacobian  matrix.  Substituting  the  expansion  into  the  governing  equations, 
the  system  becomes 


I  -  A"  — 
2 


(Qn+l-Qn)=QnAt 


(30) 


where  I  is  the  identity  matrix  and  represents  a  system  of  coupled  linear 
equations  in  the  form  of  Ax  =  b . 

Even  when  using  the  implicit  method  there  remains  a  variety  of  options 
that  finalize  the  solution  approach.  One  option  is  the  solution  of  the  equation  set 
and  here  the  LU  Decomposition  method  will  be  used.  Another  option  is  the 
formulation  of  the  Jacobian  matrix.  This  choice  itself  can  affect  the 
computational  time  required  to  perform  a  simulation.  One  approach  that  can 
greatly  reduce  run  time  is  to  determine  the  Jacobian  analytically  [10].  However, 
this  approach  requires  preprocessing  using  other  software  tools  rather  than  just 
the  tool  to  be  used  for  the  reaction  simulation.  Furthermore  this  preprocessing  is 
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required  for  each  reaction  scenario  to  be  considered.  Therefore  for  the  current 
study  this  approach  will  not  be  used. 

The  first  option  to  establishing  the  Jacobian  is  a  generic  approach  where 
the  computational  routine  uses  a  set  of  do  loops  with  the  parameters  being  set 
by  an  input  file.  Such  an  approach  does  not  require  extensive  preprocessing 
before  executing  a  simulation  but  only  reading  a  simple  file  that  defines  the 
mechanism  such  as  Figures  96,  97,  and  98.  This  approach  will  be  one  of  the 
methods  evaluated  here. 

Algorithms  utilizing  do  loops  can  be  computationally  expensive  because  of 
the  many  looping  functions.  To  overcome  this  limitation,  an  alternative  approach 
using  numerical  derivatives  will  be  evaluated.  By  calculating  the  Jacobian 
numerically,  only  a  single  looping  function  is  required  as  opposed  to  several 
embedded  functions  used  in  the  general  implicit  approach.  A  study  by  Orkwis 
and  Vanden,  which  compared  the  use  of  numerical  derivatives  and  analytic 
representation  of  the  Jacobian  terms  in  the  full  Navier-Stokes  equations,  found 
that  numerical  derivatives  can  be  used  to  accurately  represent  the  terms  of  the 
Jacobian  [31].  Assuming  that  one  has  the  function  ,  the  definition  of  its 
derivative  is 


f -/(»)- \™f(x  +  h)~f(x).  (31) 

ax  o  h 

Using  the  form  above,  the  derivatives  of  each  source  term  that  composes  the 
Jacobian  can  be  represented  as 


doj,  _  coi [paj  +  h)+  {paj )  (32) 

dpctj  h 

where  i  and  j  represent  each  species  in  the  mechanism  and  h  a  deviation  in  the 
concentration  of  species  j.  Because  the  source  term  at  state  n  has  to  be 
calculated  for  any  method  used  to  solve  the  reaction  equations,  the  numerical 
derivative  approach  requires  only  an  additional  calculation  which  can  be  done 
within  a  single  looping  function. 

The  numerical  derivatives  contain  two  sources  of  error,  truncation  error 
found  in  any  difference  equation  and  roundoff  error.  The  roundoff  error  is 
inherent  in  any  computer  model  due  to  machine  accuracy.  If  the  size  of  h  is 
outside  the  precision  of  the  computer  or  programming  environment,  the 
machine's  floating  point  representation  may  not  accurately  represent  the  effect 
of  h.  The  worst  case  of  this  is  where  the  value  of  h  is  lost  in  rounding  and  the 
resulting  derivative  becomes  zero  or  undefined. 
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6.2  Results 

As  a  starting  point  for  the  analysis,  initial  conditions  were  chosen  to  be 
representative  of  temperatures  and  pressures  that  would  be  found  in  a  high 
speed  reacting  flow  such  as  a  high  explosive  process.  For  this,  a  temperature  of 
1600  K  and  a  pressure  of  1125  kPa  will  be  used  as  a  baseline  for  comparison. 
Each  scheme  was  run  using  a  constant  time  step  of  IE-9  seconds  and  run  for  the 
three  reaction  mechanisms.  This  small  value  was  chosen  to  ensure  stability  for  all 
methods. 

Figures  99,  100,  and  101  show  plots  of  temperature  versus  time  for  each 
method  and  the  various  reaction  mechanism.  The  temperature  histories 
produced  by  the  general  explicit,  Runge-Kutta,  and  numerical  derivative  schemes 
are  almost  identical  but  the  history  of  the  general  implicit  scheme  differs. 
Although  the  general  implicit  method  deviates,  it  reaches  the  same  steady  state 
conditions  as  the  explicit  methods  and  the  numerical  derivative  method.  The 
deviation  of  the  generic  implicit  method  can  be  attributed  to  roundoff  error 
accumulated  in  the  calculations  which  has  been  found  in  other  studies  [31].  The 
19-step  and  32-step  mechanisms  produce  similar  results  but  as  the  mechanism 
becomes  more  complex  the  difference  between  the  general  explicit  results  and 
the  others  becomes  smaller. 

Another  obvious  result  is  the  difference  in  the  predicted  heat  release  as 
the  difference  mechanisms  are  chosen.  The  19-step  and  32-step  mechanisms 
results  are  somewhat  similar  but  the  7-step  is  considerably  different.  This 
reinforces  the  point  that  when  a  reduced  mechanism  is  used  fundamentally 
different  predictions  can  result.  This  is  particularly  true  if  the  reduction  was 
meant  for  a  particular  operational  scenario  but  the  problem  of  interest  falls 
outside  this  scenario.  It  further  emphasizes  the  need  to  have  better  numerical 
techniques  that  allow  for  more  complex  mechanisms  to  be  used,  particularly 
when  addressing  problems  where  the  chemical  processes  are  not  well  defined. 


Figure  99.  Temperature  history  for  each  numerical  scheme  using  7-step 
mechanism  with  a  constant  time  step  (At  =  le-9). 


Figure  100.  Temperature  history  for  each  numerical  scheme  using  19 
step  mechanism  with  a  constant  time  step  (At  =  \e-9) . 
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Figure  101.  Temperature  history  for  each  numerical  scheme  using  32- 
step  mechanism  with  a  constant  time  step  (a/  =  ie-9). 

In  this  first  set  of  results,  the  governing  equations  were  integrated  in  time 
using  a  constant  time  step.  This  allows  for  the  quantification  of  the 
computational  requirements  of  each  method.  Table  4  shows  the  required  run 
time  for  each  method.  All  simulations  were  conducted  on  the  same  single 
processor,  3.0GHz  desktop  computer.  The  generic  explicit  method  requiring  the 
least  run-time  for  each  mechanism  and  the  Runge-Kutta  method  also  proved  to 
be  extremely  efficient  in  the  constant  time  step.  This  is  due  to  the  fact  that  no 
matrix  operations  are  required.  However,  when  using  the  explicit  methods  small 
time  steps  must  be  taken  to  maintain  stability.  Another  possible  occurrence  is 
that  the  solution  can  produce  non-physical  solutions  such  as  having  mass 
fractions  less  than  zero.  Another  notable  result  in  Figure  102  is  the  times  for  the 
general  implicit  and  numerical  derivative  methods.  The  generic  implicit  method 
required  the  greatest  amount  of  run  time  to  model  the  reactions  due  to  the 
number  of  looping  functions  required  to  calculate  the  Jacobian  terms.  Results 
also  show  that  the  time  saved  by  evaluating  the  Jacobian  terms  numerically 
grows  as  the  mechanism  becomes  more  complex.  The  increase  between  the  two 
methods  grew  from  457%  for  the  7-Step  mechanism  to  1383%  for  the  32-Step 
mechanism.  The  numerical  derivatives  are,  on  average,  approximately  nine  times 
more  efficient  than  the  generic  implicit  solution  when  run  with  constant  time 
steps. 
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Reaction  Mechanism 

Scheme 

7 -step 

1 9- step 

3  2- step 

General  Explicit 

335.00 

546.00 

821.00 

Runge-Kutta 

625.00 

1555.00 

3033.00 

Numerical  Derivatives 

687.00 

2118.00 

5982.00 

General  Implicit 

3845.00 

19272.00 

82715.00 

Figure  102.  Computational  times  in  seconds  using  (At  =  \e-9)  for  each 
reaction  mechanism  and  solution  scheme. 

The  true  improvement  and  efficiency  of  the  different  schemes  comes  from 
the  use  of  variable  time  steps.  Unless  the  time  step  in  the  simulation  is  limited  to 
resolve  a  process  of  interest,  the  desire  is  to  integrate  the  governing  equations 
through  the  time  period  as  rapidly  as  possible.  Therefore  in  the  following 
simulations  each  scheme  was  allowed  to  set  its  At  based  on  stability 
requirements  and  to  keep  the  solution  out  of  non-physical  regions.  A  maximum 
allowable  value  for  At  was  set  to  be  IE-5  seconds.  Each  reaction  mechanism 
was  again  modeled  with  initial  conditions  of  T=1600K  and  P=1125kPa  and  out  to 
a  time  of  one  millisecond. 

Figures  103,  104,  and  105  show  temperature  histories  for  each  scheme 
and  mechanism  using  a  variable  time  step.  In  order  to  retain  accuracy  the 
explicit  methods  had  to  be  forced  to  take  smaller  time  steps  at  the  beginning  of 
the  calculations.  The  generic  explicit  results  using  a  variable  time  step  matches 
the  results  when  the  constant  time  step  is  employed.  The  numerical  derivatives 
and  generic  implicit  method  have  shapes  similar  to  that  of  their  constant  time 
step  solutions  with  the  numerical  derivative  approach  again  more  closely  tracking 
the  explicit  results.  This  is  due  in  part  to  the  fact  that  the  implicit  method  using 
the  numerical  derivatives  was  found  to  take  smaller  time  steps  initially.  The 
generic  implicit  time  history  contains  fewer  points,  thus  it  gives  a  rougher 
solution  than  the  numerical  derivatives  method.  This  could  easily  be  corrected  by 
further  restricting  the  time  step  of  the  generic  implicit  scheme  to  better  capture 
the  profile. 
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Figure  103.  Temperature  history  for  each  numerical  scheme  using  7- 
step  mechanism  with  a  variable  time  step. 


Figure  104.  Temperature  history  for  each  numerical  scheme  using  19- 
step  mechanism  with  a  variable  time  step. 
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Figure  105.  Temperature  history  for  each  numerical  scheme  using  32- 
step  mechanism  with  a  variable  time  step. 


Reaction  Mechanism 

Scheme 

7 -step 

1 9- step 

3  2- step 

General  Explicit 

77.00 

163.00 

571.00 

Runge-Kutta 

1464.00 

2443.00 

4987.00 

Numerical  Derivatives 

3.000 

4.000 

5.000 

General  Implicit 

4.000 

6.000 

12.000 

Figure  106.  Computational  times  in  seconds  using  a  variable  time  step 
for  each  reaction  mechanism  and  solution  scheme. 

The  run  time  results  for  each  method  when  using  a  variable  time  step  to 
solve  the  various  mechanisms  are  shown  in  Figure  106.  The  overall  run  times  are 
reduced  due  to  the  fact  that  the  schemes  are  allowed  to  take  the  maximum  time 
step  they  can  and  maintain  a  stable  and  physically  possible  solution.  Employing 
the  numerical  derivatives  is  found  to  reduce  the  computational  time  and  the 
advantage  improves  as  the  mechanism  complexity  increases.  The  advantage  of 
using  an  implicit  method  is  clear  and  can  be  attributed  directly  to  the  increased 
stability  allowing  for  larger  integration  time  steps.  It  should  be  noted  that  some 
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of  the  advantage  will  be  diminished  if  the  maximum  allowable  time  step  were 
reduced  out  of  necessity  to  resolve  important  intermediate  stages  in  the  reaction 
process. 

At  the  center  of  simulating  reacting  flow  scenarios  is  the  solution  of  the 
equations  governing  the  combustion  process.  The  two  primary  classes  of 
methods  are  explicit  and  implicit.  The  explicit  methods  require  much  less 
computational  resources  per  time  integration  step.  However,  due  to  stability 
requirements  the  time  step  when  these  methods  are  used  have  to  be  limited  to 
small  values.  This  is  acceptable  if  a  small  step  is  desired  to  ensure  resolution  of 
intermediate  chemical  processes.  However,  it  is  more  desirable  to  have  an 
integration  method  that  is  not  limited  primarily  by  the  numerical  aspects. 

Implicit  methods  offer  advantages  in  terms  on  their  stability  and  due  to 
their  nature  the  solutions  tend  to  allow  for  larger  time  integration  steps.  Though 
the  implicit  schemes  require  considerably  more  computational  time  per  step, 
these  advantages  overcome  the  shortfalls  when  considering  the  total  time  of  the 
simulation.  The  typical  implicit  method  uses  several  looping  routines  to  construct 
the  Jacobian  matrix  needed  for  the  operations.  By  introducing  numerical 
derivatives  the  resulting  method  is  found  to  produce  shorter  run  times  and 
reduce  the  roundoff  error.  These  initial  results  support  further  investigation  into 
using  numerical  derivatives  in  the  solution  of  the  reaction  equations  found  in 
combustion  flow  problems. 

7.  Results 

Initially  this  project  was  utilizing  a  code  titled  cutter  provided  by  the  US 
Air  Force  Research  Lab's  Munition  Directorate.  We  eventually  had  to  shift  to 
using  an  alternative  computational  framework  to  be  able  to  pursue  some  of  the 
topics  of  interest.  Using  the  new  framework  to  include  the  adaptive  meshing  the 
various  detonation  modeling  strategies  were  evaluated.  The  case  of  an  impact 
sending  a  shock  into  the  explosive  material  was  the  scenario  used  and  is 
depicted  in  Figure  107.  This  is  patterned  after  the  flyer  plate  problem  that  has 
been  used  to  test  the  ignition-growth  model  in  cutter  [  12].  Here  a  slug  of  air  is 
used  to  impact  the  explosive  since  no  constitutive  models  for  inerts  have  been 
integrated  into  the  current  modeling  framework. 
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\\  =  5.6  kni/s,  2.8  km/s 


Figure  107.  Configuration  used  in  the  shock  initiation  simulations. 

To  simulate  these  cases,  the  JWL  equation  of  state  was  used  for  the 
unreacted  explosive  material  and  the  detonation  products  with  the  parameters  in 
Table  1.  The  first  simulation  cases  used  the  ignition  and  growth  model  by  [23] 
and  demonstrated  by  [12].  The  parameters  for  the  reaction  are  listed  in  Table  2. 
Two  impact  cases  were  run  with  impact  velocities  of  2.8  km/s  and  5.6  km/s. 
These  two  were  selected  because  the  lower  velocity  is  known  to  not  produce  a 
self  sustaining  detonation  but  the  higher  velocity  will  [11]. 


Un  reacted 
Explosive 

Detonation 

Products 

po  (g/cc) 

1.624 

- 

A  (GPa) 

17101 

673.1 

B  (GPa) 

-3.745 

21.988 

9.8 

5.4 

r2 

.98 

1.8 

CO 

.5675 

.3 

cv 

(GPa/K) 

2.70386e-3 

1.0e-3 

E0  (GPa) 

0 

7.0 

Table  1.  J  WL  EOS  Parameters  for  the  TNT  and  detonation  products. 
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Reaction 

Rate 

Parameters 

I=5.0e7 

a=0 

b=.667 

O 

'xt1 

II 

X 

Gj=3.6e8 

y=1.2 

O 

i — 1 

II 

u 

d=.667 

G2=  0 

e=1.0 

g=.lll 

O 

i — 1 

II 

IN 

F  =0.03 

igmax 

O 

i — 1 

II 

3 

s 

k5 

^Glmin  ^-0 

Table  2.  J  WL  EOS  Parameters  for  the  TNT  and  detonation  products. 

7.1  Single  Step  Models 

The  first  reaction  models  tested  were  those  that  use  a  single  step  to 
represent  the  reactions  involved  in  the  detonation.  Figures  108,  109,  and  110 
are  a  sequence  of  images  showing  the  predicted  response  of  the  TNT  using  the 
ignition  and  growth  when  impacted  at  2.8  km/s.  A  compression  wave  is  seen  to 
travel  through  the  energetic  material  and  associated  heating  occurs.  However,  a 
self-sustaining  detonation  is  not  produced.  For  this  simulation  only  one  level  of 
cells  was  used.  Figures  111,  112,  and  113  show  the  response  when  the  impact 
velocity  is  raised  to  5.6  km/s.  In  the  pressure  data,  the  von  Neumann  pressure 
behind  the  shock  of  the  detonation  front  is  evident  and  the  value  is  consistent 
with  the  experimental  data.  The  detonation  forms  and  is  self  sustained  through 
the  full  length  of  the  explosive. 

Figures  114,  115,  and  116  show  the  5.6  km/s  case  results  when  two  levels 
of  computational  cells  are  used.  The  parameter  e  plotted  with  the  velocity  is  the 
adaptation  parameter  (z)  discussed  earlier.  Using  two  levels  does  not  offer 
much  advantage  for  the  ID  impact  problem  but  was  tested  here. 

The  results  using  the  ignition  and  growth  model  demonstrate  correct 
response  to  impact  and  produce  the  correct  data  in  terms  of  post  detonation 
pressure.  However,  the  temperatures  determined  using  the  JWL  EOS  data  are 
not  correct.  The  detonation  temperature  for  TNT  is  known  to  be  approximately 
3700  K.  This  is  not  a  surprise  since  the  JWL  EOS  was  not  formulated  with  the 
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intention  to  correctly  predict  temperature.  The  ignition  and  growth  model  does 
not  even  use  this  variable  in  its  calculation.  However,  this  is  a  shortcoming  if  a 
single  reaction  model  is  desired  to  predict  all  possible  reaction  scenarios. 


Figure  108.  Sequence  of  images  showing  the  results  using  the  ignition 
and  growth  model  and  a  single  level  of  cells  for  the  K  =  2Mm/s  case 

shortly  after  impact. 
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Figure  109.  I  mage  showing  the  results  using  the  ignition  and  growth 
model  and  a  single  level  of  cells  for  the  f  '  =  2.8 km/s  case  some  time  after 

impact.  No  self  sustaining  detonation  has  formed. 
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Figure  110.  I  mage  showing  the  results  using  the  ignition  and  growth 
model  and  a  single  level  of  cells  for  the  K  =  2.8 km/s  case  well  after 

impact.  No  self  sustaining  detonation  has  formed. 


Figure  111.  Sequence  of  images  showing  the  results  using  the  ignition 
and  growth  model  and  a  single  level  of  cells  for  the  l';  =  5 .6km/s  case 

shortly  after  impact. 


47 


Figure  112.  I  mage  showing  the  results  using  the  ignition  and  growth 
model  and  a  single  level  of  cells  for  the  (  ’  =  5.6km/s  case  some  time  after 

impact.  A  self  sustaining  detonation  has  formed. 


Figure  113.  I  mage  showing  the  results  using  the  ignition  and  growth 
model  and  a  single  level  of  cells  for  the  (’  =  5 ,6km/ s  case  well  after 

impact  and  the  formation  of  the  self  sustaining  detonation. 
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Figure  114.  Sequence  of  images  showing  the  results  using  the  ignition 
and  growth  model  and  two  levels  of  cells  for  the  Vi  =  5.6km/s  case 

shortly  after  impact. 


Figure  115.  I  mage  showing  the  results  using  the  ignition  and  growth 
model  and  two  levels  of  cells  for  the  K  =  5 Mm/s  case  some  time  after 

impact.  A  self  sustaining  detonation  has  formed. 
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Figure  116.  I  mage  showing  the  results  using  the  ignition  and  growth 
model  and  two  levels  of  cells  for  the  Vi  =  5.6 kmls  case  well  after  impact 

and  the  formation  of  the  self  sustaining  detonation. 


Figure  117. 1  mages  showing  the  results  using  the  hot  spot  and  balance 
model  and  a  single  level  of  cells  for  the  (’  =  2.8 kmls  case  shortly  after 

impact. 
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Figures  117  and  118  are  a  sequence  of  images  showing  the  predicted 
response  of  the  TNT  using  the  hot  spot  and  balance  model  when  impacted  at  2.8 
km/s.  Figures  119,  120,  and  121  show  the  results  for  the  5.6  km/s  case.  The 
results  are  seen  to  be  similar  to  those  from  the  ignition  and  growth  model.  Again 
the  pressure  and  wave  dynamics  are  correctly  captured  but  the  temperature 
profiles  are  not  accurate.  The  hot  spot  and  balance  model  is  a  temperature 
based  reaction  model  however  since  the  material  model  does  not  accurately 
predict  temperature  in  the  shocked  material,  the  model  is  not  that  sensitive  to 
the  actual  temperature.  This  insensitivity  is  clear  if  the  formation  of  the  model  is 
reviewed.  The  temperature  dependence  is  derived  by  using  the  known 
detonation  pressure  and  induction  time  for  the  material.  In  the  work  by  [17]  the 
Mie-Gruneisen  equation  of  state  was  used  so  the  calculated  temperature  is 
different  and  is  reported  to  be  more  consistent  with  the  know  detonation 
temperature. 


Figure  118. 1  mages  showing  the  results  using  the  hot  spot  and  balance 
model  and  a  single  level  of  cells  for  the  Vi  =  2.8 km/s  case  some  time  after 

impact. 
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Figure  119. 1  mages  showing  the  results  using  the  hot  spot  and  balance 
model  and  a  single  level  of  cells  for  the  Vt  =  5 ,6kmls  case  shortly  after 

impact. 
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Figure  120. 1  mages  showing  the  results  using  the  hot  spot  and  balance 
model  and  a  single  level  of  cells  for  the  f  '  =  5.6km/s  case  some  time  after 


impact.  A  self  sustaining  detonation  has  formed. 
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Figure  121. 1  mages  showing  the  results  using  the  hot  spot  and  balance 
model  and  a  single  level  of  cells  for  the  (’  =  5 ,6km/s  case  well  after 

impact  and  the  formation  of  the  self  sustaining  detonation. 
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7.2  Multi-step,  Finite- rate  Model 

The  single-step  kinetics  models  perform  adequately  when  the  interest  is 
to  reproduce  the  travel  of  the  detonation  wave.  However,  they  do  not  offer  a 
wide  range  of  applicability  and  will  not  accurately  model  all  scenarios  of  interest. 
As  mentioned,  one  objective  of  this  project  was  to  determine  the  feasibility  of 
having  a  kinetics  based  model  that  could  correctly  represent  a  wide  range  of 
scenarios  involving  energetic  materials.  The  range  of  possible  scenarios  run  from 
slow  events  such  as  cook-off  to  the  case  of  detonations.  The  recent,  successful 
work  in  using  multi-step,  finite-rate  models  in  simulating  cook-off  [40]  provided  a 
good  point  of  departure. 

In  this  earlier  study,  the  decomposition  of  the  explosive  RDX  was 
represented  using  a  three-step,  four-species  model  of  the  form 


A  — »  B,rt  =  Ztexp(-  EJRT )pA 

(33) 

B  — »  C,r2  =  Z2exp(-  E2/RT )pB 

(34) 

C  — ^  D,  v2  =  Z3exp(-  EJRT  )p1c 

(35) 

where  A  and  B  are  solid  species,  and  C  and  D  are  gases.  The  parameter  r  is  the 
reaction  rate  and  z,  is  the  rate  factor  and  Ej  is  the  activation  energy  for 

reaction  i.  This  form  is  similar  to  the  Arrenhius  rate  law  typically  used  in  gas 
phase  combustion  modeling.  The  parameter  p]  is  the  mass  fraction  of  species 

j .  The  parameters  for  the  RDX  decomposition  can  be  found  in  Table  3.  Here  this 
baseline  model  will  be  used  for  the  detonation  case  to  see  how  it  performs.  Prior 
to  that,  the  one-dimensional-time-to-explosion  (ODTX)  case  for  the  RDX 
simulated  by  [40]  was  recreated. 


Reaction 

step 

In  zt 

Ek  (kj/g 
mole  K) 

qk  (J/g) 

A^B 

43.70s-1 

197.2 

139.0 

B^C 

38.90s-1 

184.7 

-535.6 

C^D 

32.69 cnv’lsg 

142.8 

-4680 

Table  3.  Reaction  kinetics  parameters  for  the  decomposition  of  RDX. 

The  ODTX  scenario  is  modeled  using  the  chemical  reaction  source  model 
discussed  earlier  and  the  time  to  explosion  was  determined  for  a  range  of 
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temperatures.  In  the  simulations  the  fluid  dynamics  operations  can  be  ignored 
and  a  mass  of  explosive  is  suddenly  increased  to  the  specified  temperature.  The 
time  history  of  the  species  concentrations  and  the  temperatures  were  recorded 
and  presented  in  Figure  122.  Figure  123  summarizes  the  predicted  times  to 
explosion  as  a  function  of  temperature.  Also  shown  in  the  figure  are  the  results 
from  the  earlier  study  [40]  and  experimental  data.  The  current  results  show 
shorter  times  to  explosion  for  the  higher  temperatures  than  those  predicted  in 
the  earlier  study  and  the  few  experimental  points  in  that  region.  However,  at  the 
lower  temperatures  the  current  results  compare  well  with  the  experimental  data. 


the  one-dimensional-time-to-explosion  (ODTX)  simulations. 
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Figure  123.  Calculated  ODTX  times  using  the  current  model  compared 
to  experimental  data  and  earlier  studies. 

This  three-step,  four-species  model  was  tested  within  the  context  of  the 
impact  scenario  run  with  the  one-step  models.  The  5.6  km/s  case  was  run  and 
sample  results  are  shown  in  Figures  124,  125,  and  126.  The  solution  shows  that 
a  shock  is  driven  through  the  energetic  material  by  the  reaction  front,  however 
the  pressure  does  not  reach  the  von  Neuman  value  and  the  front  speed  is  not 
equal  to  the  known  detonation  velocity.  This  is  not  a  surprise  since  the  reaction 
paths  and  time  scales  for  the  detonation  regime  will  be  different  for  the  cook-off 
regime  for  which  the  current  model  was  developed.  In  the  next  section  some 
initial  work  will  be  presented  related  to  modifying  the  mechanism  to  also  cover 
the  detonation  regime. 
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Figure  124.  Images  showing  the  results  using  the  multi-step,  finite- 
rate  model  and  a  single  level  of  cells  for  the  Vi  =  5 .6km/ s  case  just  after 

impact. 


Figure  125.  Images  showing  the  results  using  the  multi-step,  finite- 
rate  model  and  a  single  level  of  cells  for  the  K  =  5.6km/s  case  shortly 

after  impact. 
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rate  model  and  a  single  level  of  cells  for  the  /  '  =  5.6km/s  case  some  time 
after  impact. 

In  the  earlier  work  by  [40],  a  one-step  model  was  used  to  simulate  the 
travel  of  the  detonation  after  initiation  from  the  thermal  explosion  process.  The 
model  was  used  here  to  simulate  the  impact  scenario  for  which  the  other  one- 
step  models  were  used.  Since  this  model  is  designed  for  the  detonation  regime 
just  the  5.6  km/s  case  was  simulated  and  the  results  are  shown  in  Figures  127 
and  128.  The  model  is  seen  to  produce  the  correct  von  Neumann  pressure. 
Figure  130  compares  the  predicted  time  histories  of  pressure  at  a  location  in  the 
explosive  for  the  detonation  case  for  each  one-step  model  and  the  three-step 
model.  The  one-step  models  provide  similar  predictions  of  the  pressure  profiles 
but  the  arrival  times  vary. 

The  three-step  model  is  found  to  not  produce  the  von  Neumann  spike 
known  to  exist.  This  may  be  due  to  the  fact  that  since  the  reaction  process  in  the 
detonation  is  being  represented  by  a  multi-step  model,  the  required  resolution  is 
extremely  rigid.  This  can  be  overcome  using  the  adaptation  methods  detailed 
here.  However,  it  would  be  desirable  to  have  a  model  that  might  require  grid 
adaptation  for  exact  resolution  but  not  be  totally  dependent  on  the  grid  to  show 
the  general  nature  of  the  phenomena.  It  is  clear  that  all  1-step  models  do 
reproduce  the  general  trend.  A  strategy  to  combine  the  characteristics  of  the 
multi  and  single  step  models  will  be  presented  in  the  next  section. 
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Figure  127. 1  mages  showing  the  results  using  the  hot  spot  and  balance 
model  and  a  single  level  of  cells  for  the  !'  =  5.6 km/s  case  shortly  after 

impact. 


Figure  128. 1  mages  showing  the  results  using  the  hot  spot  and  balance 
model  and  a  single  level  of  cells  for  the  l  '  =  5.6km/s  case  some  time  after 

impact.  A  self  sustaining  detonation  has  formed. 
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Figure  129. 1  mages  showing  the  results  using  the  hot  spot  and  balance 
model  and  a  single  level  of  cells  for  the  l)  =  5 ,6km/s  case  well  after 

impact  and  the  formation  of  the  self  sustaining  detonation. 


Figure  130.  Comparison  of  the  predicted  pressure  profiles  for  the 
detonation  case  using  the  various  one-step  models. 
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7.3  Modified  Multi-step,  Finite- rate  Model 

As  shown  here,  a  single-step  reaction  model  is  sufficient  to  represent  a 
detonation  process.  Though  the  rate  based  models  are  dependent  on  local 
conditions  a  multi-step  model  is  needed  to  represent  the  full  spectrum  of 
reaction  scenarios,  an  objective  of  this  project.  As  shown,  the  multi-step  models 
like  that  used  by  [40]  can  represent  these  other  processes.  However,  it  is  found 
to  not  easily  show  the  detonation  case  and  when  used  in  the  earlier  study  a 
switch  is  made  from  the  multi-step  model  to  a  simpler  1-step  model  once  an 
explosion  occurs.  This  switching  process  would  be  difficult  to  make  autonomous 
and  an  elaborate  method  would  be  needed  if  the  detonation  mode  needed  to  be 
switched  back  to  the  deflagration  case.  A  universally  application  approach  may 
be  possible  by  incorporating  characteristics  used  in  other  combustion  problems. 

Multi-step,  finite-rate  reaction  models  are  routinely  used  to  model  a  wide 
range  of  combustion  processes.  With  enough  steps,  multiple  reaction  paths  are 
adequately  represented.  By  representing  the  multiple  reaction  paths,  a  full  range 
of  combustion  scenarios  can  be  simulated.  One  example  is  how  mechanisms  can 
represent  the  explosion  limits  of  fuel-oxidizer  mixtures.  An  example  can  be  seen 
in  the  hydrogen-air  mechanisms  in  Figures  96  and  98.  Though  there  is  a 
difference  of  25  steps,  both  reproduce  slow  burning  and  explosion  scenarios. 
This  can  be  shown  by  looking  at  results  from  an  earlier  study  ([10])  using  the 
two  mechanisms  in  simulating  shock  initiated  combustion,  similar  to  a 
detonation. 

Figure  131  shows  an  experimental  image  of  the  scenario  which  is  a 
projectile  shot  into  a  mixture  of  hydrogen  and  air.  The  curved  shock  is  visible. 
Also  shown  is  the  computed  results  of  induction  time  of  the  mixture  depending 
on  which  part  of  the  shock  the  mixture  crosses.  The  conditions  on  the  respective 
streamline  behind  the  shock  (point  2  in  the  insert  in  Figure  131)  are  used  to 
integrate  the  reaction  equations.  At  locations  of  the  curved  shock  when  the 
shock  angle  is  greater  than  approximately  50°,  relative  to  the  streamline,  the 
induction  time  is  long.  Otherwise  they  are  very  short. 
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Figure  131.  Shock  initiated  combustion  scenarios  and  predicted 
induction  times. 

In  the  hydrogen-oxygen  example,  the  explosion  limit  is  controlled  by  the 
branching  reaction 


H  +  02^0H  +  0 


(36) 


and  the  terminating  reaction 


H  +  02+M  <^H02+M  (37) 

where  M  represents  any  compound  in  the  mixture.  These  third  body  reactions 
are  highly  dependent  on  pressure.  Therefore,  the  local  conditions  to  include 
pressure  dictate  whether  a  slow  combustion  process  occurs  or  a  detonation  is 
produced  by  the  model.  A  similar  approach  is  possible  for  use  for  the  current 
energetics  problems.  To  test  this  idea,  the  multi-step  model  represented  by 
Equation  (33),  (34),  and  (35)  is  augmented  with  the  additional  step 

A^D,rA=  Z4exp(-  EJRT  )pA  (38) 

which  represents  an  alternative  reaction  path.  Figure  132  shows  the  relative 
rates  at  which  each  reaction  occurs  as  temperature  increases.  The  rate  of  the 
step  represented  in  Equation  (38)  was  set  such  that  it  would  be  dominate  at  the 
higher  temperatures.  The  reaction  step  could  also  be  fashioned  to  be  dependent 
of  pressure. 
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Figure  132.  Comparison  of  the  rates  for  each  step  in  the  modified 
model. 

This  4-step  model  is  compared  to  the  original  3-step  model  by  seeing  how 
it  responds  to  the  same  initial  conditions.  This  is  done  by  comparing  the 
predicted  time  histories.  One  way  is  to  plot  the  mass  fraction  of  each  material  as 
a  function  of  time  as  shown  in  Figure  133.  Another  useful  comparison  can  be 
made  by  plotting  various  material  concentrations  versus  others.  In  Figure  133 
the  intermediate  products  B  and  C  are  plotted  as  a  function  of  the  final  product 
D.  The  3-step  and  4-step  results,  which  are  identical,  of  course  start  with  D=0 
and  terminate  when  D=l.  The  symbols  denote  intervals  in  the  integration  and 
are  evenly  spaced  in  time.  The  reaction  is  seen  to  process  slowly  with  B  and  C 
increasing  up  to  a  point  and  it  then  rapidly  converges  to  the  final  state. 

Similar  results  are  seen  when  the  the  initial  temperature  is  set  to  600K  as 
shown  in  Figure  134.  However,  the  slope  early  in  the  process  and  the  location  of 
the  transition  changes.  But  when  the  initial  temperature  is  set  to  1000K,  there  is 
a  distinct  difference  in  the  3-step  and  4-step  results.  As  evident  in  Figure  135, 
the  3-step  model  shows  results  similar  to  the  lower  temperatures  but  the 
transition  occurs  much  earlier.  The  highly  nonlinear  nature  of  the  reaction 
process  is  one  reason  the  integration  of  the  multi-step  models  are  much  more 
difficult.  When  the  4th  step  in  added  to  the  model,  this  higher  initial  temperature 
causes  the  reaction  to  immediately  move  through  the  transition  state  to  the  final 
state.  This  is  due  solely  to  the  alternate  reaction  path  introduced  to  the  model. 
Figure  136  shows  a  variety  of  comparisons  by  plotting  various  material 


63 


concentrations  versus  others.  Introducing  these  alternative  reaction  paths  appear 
to  be  a  valid  strategy  for  having  a  single  reaction  model  that  covers  all  possible 
reaction  scenarios. 


Figure  133.  Comparison  of  the  3-step  and  4-step  reaction  model  results 
with  an  initial  temperature  of  500K. 


Figure  134.  Comparison  of  the  3-step  and  4-step  reaction  model  results 
with  an  initial  temperature  of  600K. 


Figure  135.  Comparison  of  the  3-step  and  4-step  reaction  model  results 
with  an  initial  temperature  of  1000K.  The  plot  showing  concentration 
as  a  function  of  time  is  from  the  baseline  3-step  model  since  the  4-step 
transitions  immediately  to  the  final  state. 
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Figure  136.  A  variety  of  species  space  plots  from  the  3-step  and  4-step 
reaction  model  results  with  an  initial  temperature  of  1000K. 


7.4  Example  Full  Weapon  System  Simulation 

The  ultimate  objective  is  to  have  a  modeling  strategy  that  correctly 
simulates  a  variety  of  munition  problems.  One  example  is  the  case  of  premature 
burning  of  energetic  materials  in  a  weapon  system  due  to  mechanical  stresses 
during  employment.  Such  a  situation  could  cause  the  weapon  to  not  provide  the 
full  impact  desired  when  initiation  occurs  since  some  of  the  energy  would  have 
already  been  released  at  a  rate  slower  than  the  detonation  case.  The  methods 
proposed  here  can  be  coupled  together  to  provide  such  a  modeling  tool. 

To  demonstrate  the  use  of  the  multi-step  reaction  model  with  the 
adaptive  meshing  framework,  a  weapon  system  scenario  was  simulated.  The 
system  is  a  typical  penetrating  warhead  filled  with  explosive  material.  To 
demonstrate  the  strategy  a  portion  of  the  explosive  bed  was  heated  at  a 
temperature  that  would  be  produced  during  a  mechanical  stressing  process. 
Figure  137  shows  a  sequence  of  images  showing  the  temperature  and  pressure 
profile  along  a  plane  through  the  center  of  the  item.  The  general  process  is 
captured  and  the  progression  of  the  burning  is  tracked.  The  simulation  was  not 
run  to  the  point  of  detonation.  To  fully  simulate  this  scenario  an  additional 
thermal  model  to  include  heat  transfer  through  the  energetic  bed  would  be 
needed. 
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Figure  137.  Sequence  of  images  showing  the  solution  of  the  three-step 
reaction  model  for  a  full  system  simulation. 


8.  Findings  and  Future  Work 

Kinetics  based  modeling  of  detonations  is  a  viable  and  reliable  simulation 
option.  Here  a  variety  of  1-step  model  types  have  been  investigated  and  all 
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provide  the  same  general  results.  However,  these  single  reaction  path  models  do 
not  offer  an  option  to  simulate  off  detonation  cases.  To  correctly  simulate  these 
scenarios  a  multi-step  model  is  needed. 

A  multi-step  model  previously  used  to  simulate  cook-off  type  scenarios 
was  tested  for  the  detonation  case.  It  was  found  to  not  provide  accurate  results. 
This  was  due  in  part  to  the  fact  that  the  mechanism  was  not  designed  to  cover 
this  regime.  However  it  was  shown  that  by  augmenting  the  basic  mechanism 
with  an  additional  reaction  path  the  response  to  the  detonation  case  is  more  in 
line  with  the  1-step  models.  This  allows  for  the  reaction  model  to  dynamically 
switch  depending  on  the  local  conditions.  It  also  relaxes  some  of  the  numerical 
stiffness  of  the  equation  system. 

Multi-step  reaction  models  have  been  used  for  quite  some  time  in  the 
simulation  of  gas  phase  combustion.  The  same  general  approaches  can  be 
applied  to  the  reaction  of  solid  phase  energetics.  Some  of  the  delay  in  using 
these  approaches  for  this  class  of  problems  has  been  the  numerical  aspects  of 
the  problem.  Advancing  techniques  such  as  the  numerical  derivative  approach 
shown  here  should  provide  robust  modeling  tools.  This  approach  also  allows  for 
any  reaction  mechanism  to  be  easily  integrated  into  a  simulation. 

The  multi-step  reaction  model  that  includes  various  reaction  paths 
provides  the  dynamic  model  which  was  an  initial  objective  of  this  study.  The 
models  have  been  benchmarked  against  test  data  such  as  the  ODTX  scenarios. 
An  example  model  was  shown  here  that  responds  differently  to  the  burning  and 
detonation  regime.  Future  work  should  be  performed  to  validate  a  multi-step 
model  that  incorporates  the  multiple  reaction  paths. 
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